/*
 * Class:        GammaProcess
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Pierre Tremblay, Jean-Sebastien Parent & Maxime Dion
 * @since        July 2003

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.stochprocess;
import umontreal.ssj.rng.*;
import umontreal.ssj.randvar.*;
import umontreal.ssj.probdist.GammaDist;
import umontreal.ssj.util.Num;

/**
 * This class represents a *gamma* process @cite fMAD98a&thinsp; (page 82)
 * @f$\{ S(t) = G(t; \mu, \nu) : t \geq0 \}@f$ with mean parameter
 * @f$\mu@f$ and variance parameter @f$\nu@f$. It is a continuous-time
 * process with stationary, independent gamma increments such that for any
 * @f$\Delta t > 0@f$,
 * @anchor REF_stochprocess_GammaProcess_GammaEqn
 * @f[
 *   S(t + \Delta t) = S(t) + X,\tag{GammaEqn}
 * @f]
 * where @f$X@f$ is a random variate from the gamma distribution
 * <tt>Gamma</tt>@f$(\mu^2\Delta t / \nu, \mu/ \nu)@f$.
 *
 * In this class, the gamma process is sampled sequentially using equation (
 * {@link REF_stochprocess_GammaProcess_GammaEqn GammaEqn} ).
 *
 * <div class="SSJ-bigskip"></div><div class="SSJ-bigskip"></div>
 */
public class GammaProcess extends StochasticProcess {
   // Make sure that process is strictly increasing: when two successive steps
   // are equal, reset the new one to   old*(1 + EPS)
    protected static final double EPS = 1.0e-15;
    protected boolean      usesAnti = false;

    protected RandomStream stream; // used in the generatePath(), so must be kept
    protected GammaGen     Ggen;
    protected double       mu,
                           nu;
    protected double       muOverNu,
                           mu2OverNu;

    // For precomputations for standard GP
    protected double[]     mu2dtOverNu;

   protected void setLarger (double[] path, int left, int mid, int right) {
      /* make sure that path[left] < path[mid] < path[right] by adding a
         small epsilon to the last two */
      double myEps;
      if (path[left] >= 0)
         myEps = EPS;
      else
         myEps = -EPS;
      path[mid] = path[left] * (1.0 + myEps);
      if (path[mid] >= path[right])
         path[right] = path[mid] * (1.0 + myEps);
   }


   protected double setLarger (double[] path, int left, int right) {
      /* make sure that path[left] < s < path[right] by adding a
        small epsilon to the last two; return s. */
      double myEps;
      if (path[left] >= 0)
         myEps = EPS;
      else
         myEps = -EPS;
      double s = path[left] * (1.0 + myEps);
      if (s >= path[right])
         path[right] = s * (1.0 + myEps);
      return s;
   }


   protected double setLarger (double v) {
      /* return a number that is slightly larger than v */
      if (v >= 0)
         return v*(1.0 + EPS) + Num.DBL_EPSILON;
      else
         return v*(1.0 - EPS);
    }

   /**
    * Constructs a new `GammaProcess` with parameters @f$\mu=
    * \mathtt{mu}@f$, @f$\nu= \mathtt{nu}@f$ and initial value @f$S(t_0)
    * = \mathtt{s0}@f$. The gamma variates @f$X@f$ in (
    * {@link REF_stochprocess_GammaProcess_GammaEqn
    * GammaEqn} ) are generated by inversion using `stream`.
    */
   public GammaProcess (double s0, double mu, double nu,
                        RandomStream stream) {
        this (s0, mu, nu,  new GammaGen (stream, 1.0));
// Note that the parameters (the 1.0) supplied to
// the GammaGen constructors are meaningless.
// The 'real' parameters are supplied at time of generation
    }

   /**
    * Constructs a new `GammaProcess` with parameters @f$\mu=
    * \mathtt{mu}@f$, @f$\nu= \mathtt{nu}@f$ and initial value @f$S(t_0)
    * = \mathtt{s0}@f$. The gamma variates @f$X@f$ in (
    * {@link REF_stochprocess_GammaProcess_GammaEqn
    * GammaEqn} ) are supplied by the gamma random variate generator
    * `Ggen`. Note that the parameters of the
    * @ref umontreal.ssj.randvar.GammaGen object `Ggen` are not important
    * since the implementation forces the generator to use the correct
    * parameters (as defined above).
    */
   public GammaProcess (double s0, double mu, double nu, GammaGen Ggen) {
        this.mu    = mu;
        this.nu    = nu;
        this.x0    = s0;
        this.Ggen     = Ggen;
        stream = Ggen.getStream();
    }

   public double nextObservation()  {
        double s = path[observationIndex];
        double v = s;
        s += Ggen.nextDouble (stream, mu2dtOverNu[observationIndex], muOverNu);
        if (s <= v)
             s = setLarger (v);
        observationIndex++;
        observationCounter = observationIndex;
        path[observationIndex] = s;
        return s;
    }

/**
 * Generates and returns the next observation at time @f$t_{j+1} =
 * \mathtt{nextTime}@f$, using the previous observation time @f$t_j@f$
 * defined earlier (either by this method or by
 * <tt>setObservationTimes</tt>), as well as the value of the previous
 * observation @f$X(t_j)@f$. *Warning*: This method will reset the
 * observations time @f$t_{j+1}@f$ for this process to `nextT`. The user must
 * make sure that the @f$t_{j+1}@f$ supplied is @f$\geq t_j@f$.
 */
public double nextObservation (double nextT) {
        // This method is useful for generating variance gamma processes
        double s = path[observationIndex];
        double v = s;
        double previousT = t[observationIndex];
        observationIndex++;
        observationCounter = observationIndex;
        t[observationIndex] = nextT;
        double dt = nextT - previousT;
        s += Ggen.nextDouble (stream, mu2OverNu*dt, muOverNu);
        if (s <= v)
             s = setLarger (v);
        path[observationIndex] = s;
        return s;
    }

   /**
    * Generates, returns and saves the path @f$\{X(t_0), X(t_1), …,
    * X(t_d)\}@f$. The gamma variates @f$X@f$ in (
    * {@link REF_stochprocess_GammaProcess_GammaEqn
    * GammaEqn} ) are generated using the
    * @ref umontreal.ssj.rng.RandomStream `stream` or the
    * @ref umontreal.ssj.rng.RandomStream included in the
    * @ref umontreal.ssj.randvar.GammaGen `Ggen`.
    */
   public double[] generatePath()  {
        double s = x0;
        double v;
        for (int i = 0; i < d; i++) {
            v = s;
            s += Ggen.nextDouble(stream, mu2dtOverNu[i], muOverNu);
            if (s <= v)
               s = setLarger (v);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

   /**
    * Generates, returns and saves the path @f$ \{X(t_0), X(t_1), …,
    * X(t_d)\}@f$. This method does not use the
    * @ref umontreal.ssj.rng.RandomStream `stream` nor the
    * @ref umontreal.ssj.randvar.GammaGen `Ggen`. It uses the vector of
    * uniform random numbers @f$U(0, 1)@f$ provided by the user and
    * generates the path by inversion. The vector `uniform01` must be of
    * dimension @f$d@f$.
    */
   public double[] generatePath (double[] uniform01) {
        double s = x0;
        double v;
        for (int i = 0; i < d; i++) {
            v = s;
            s += GammaDist.inverseF(mu2dtOverNu[i], muOverNu, 10, uniform01[i]);
            if (s <= v)
                s = setLarger (v);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

   /**
    * Sets the parameters @f$S(t_0) = \mathtt{s0}@f$, @f$\mu=
    * \mathtt{mu}@f$ and @f$\nu= \mathtt{nu}@f$ of the process.
    * *Warning*: This method will recompute some quantities stored
    * internally, which may be slow if called repeatedly.
    */
   public void setParams (double s0, double mu, double nu) {
        this.x0 = s0;
        this.mu = mu;
        this.nu = nu;
        if (observationTimesSet) init(); // Otherwise no need to.
}

   /**
    * Returns the value of the parameter @f$\mu@f$.
    */
   public double getMu() { return mu; }

   /**
    * Returns the value of the parameter @f$\nu@f$.
    */
   public double getNu() { return nu; }

   /**
    * Resets the  @ref umontreal.ssj.rng.RandomStream of the
    * @ref umontreal.ssj.randvar.GammaGen to `stream`.
    */
   public void setStream (RandomStream stream) {
        this.stream = stream;
        Ggen.setStream (stream);
    }

   /**
    * Returns the  @ref umontreal.ssj.rng.RandomStream `stream`.
    */
   public RandomStream getStream() { return stream; }


    protected void init() {
    // init of StochasticProcess only sets path[0] = x0 if observation times are set.
        super.init();
        muOverNu  = mu / nu;
        mu2OverNu = mu * mu / nu;
        mu2dtOverNu = new double[d];
        if (observationTimesSet) {
            for (int i = 0; i < d; i++)
                mu2dtOverNu[i] = mu2OverNu * (t[i+1] - t[i]);
        }
    }

}